c  
c  solution of Laplace's eq on rectangle
c           0 < x < xmax, 0 < y < ymax
c
c  computes data for matlab program pclape55
c  error determined using exact solution
c  
c	f95 plaplace5.f	    a.out
c  
c  pcgm used to solve linear system  (m=tridiagonal)
c  
	implicit real*8(a-h,o-z)

	open(unit=7,file="p2_iterations.txt")
	write(7,201)  
  201 format('m tridiagonal')
c  201 format(7x,'n',4x,'iterations (pcgm)',5x,'time (pcgm)')
  
	imin=3
	imax=6
	nmax=15
	ap=0.5*(imax-imin)/(nmax-1.0)
	bp=0.5*(nmax*imin-imax)/(nmax-1.0)
	
	do 80 inn=1,nmax
	pow2=ap*inn+bp
	nx=nint(10**pow2)
	ny=nx
	n=nx*ny
	
	call loop(inn,nx,ny,n,it,xmin_time)

c  
c  format statements
c  
	write(7,202)  n, it, xmin_time
  202 format(4x,i8,6x,i4,12x,e13.5)
   80 continue
	close(unit=7)
	end
	
c  
c  main loop 	
c
	subroutine loop(inn,nx,ny,n,it,xmin_time)
	implicit real*8(a-h,o-z)
	real*4 time1, time2
	parameter(xmax=1.0,ymax=1.0,tol=0.000001)
	dimension b(n),v(n),q(n)
	dimension r(n),d(n),z(n),sol(n)

	dx=xmax/(nx+1)
	dy=ymax/(ny+1)
	bb=dy*dy/(dx*dx)
	
	call fourier(n,nx,ny,dx,dy,sol)
	call mult(n,nx,ny,bb,sol,b)
	
	iloops=10
	if(n.gt.1e5)  iloops=5
	
	do 70 in=1,iloops
	
	call cpu_time( time1 )
c  
c  define starting value for v
c  
	v=0.0
	call mult(n,nx,ny,bb,v,q)
	r=b-q
	call solve2(n,nx,ny,dx,dy,bb,r,z)
	d=z
	rz=dot_product(z,r)
c  
c  pcgm iteration
c  
	do 20 it=1,5000
		call mult(n,nx,ny,bb,d,q)
		dad=dot_product(d,q)
		alpha=rz/dAd
		v=v+alpha*d
		r=r-alpha*q
		call solve2(n,nx,ny,dx,dy,bb,r,z)
		rz0=rz
		rz=dot_product(z,r)
		err1=maxval(abs(sol-v))
		if(err1.lt.tol) go to 22	
		beta=rz/rz0
		d=z+beta*d
   20   continue
   22	call cpu_time( time2 )
	elapsed_time=time2-time1
	if(in.eq.1)  xmin_time=elapsed_time
	if(in.ne.1.and.elapsed_time.lt.xmin_time) xmin_time=elapsed_time
   70 continue
	return
	end
c  
c  calculate q = Ap
c  
	subroutine mult(n,nx,ny,bb,p,q)
	implicit real*8(a-h,o-z)
	dimension q(n),p(n)
	np=n-nx
	q=2.0*(1.0+bb)*p
	do 10 i=1,n
        	if(i.gt.1) q(i)=q(i)-bb*p(i-1)
        	if(i.lt.n) q(i)=q(i)-bb*p(i+1)
        	if(i.le.np) q(i)=q(i)-p(i+nx)
        	if(i.gt.nx) q(i)=q(i)-p(i-nx)
   10   continue
	do 20 j=1,ny-1
		i=j*nx
		q(i)=q(i)+bb*p(i+1)
		q(i+1)=q(i+1)+bb*p(i)
   20	continue	
	return
	end
c  
c  solve Mz=r
c  
	subroutine solve2(n,nx,ny,dx,dy,bb,r,z)
	implicit real*8(a-h,o-z)
	dimension r(n),v(n),z(n)
	dd = 2.0*(1.0+bb)
	w = dd
	z(1) = r(1)/w
	v(1)=-bb/w
	do 10 i=2,n
		w = dd + bb*v(i-1)
    		v(i) = - bb/w
    		z(i) = ( r(i) + bb*z(i-1) )/w
   10	continue
	do 20 j=n-1,1,-1
   		z(j) = z(j) - v(j)*z(j+1)
   20	continue
	return
	end
c  
c  specify boundary conditions
c      
	function gT(x)
	implicit real*8(a-h,o-z)
	gT=0.0
	if(x.ge.0.25.and.x.le.0.75) gT=1.0
	return
	end
      
      function gB(x)
      implicit real*8(a-h,o-z)
      gB=0
      return
      end
      
      function gR(y)
      implicit real*8(a-h,o-z)
      gR=0.0
      return
      end
      
      function gL(y)
      implicit real*8(a-h,o-z)
      gL=0.0
      return
      end
c  
c  calculate exact solution
c
	subroutine fourier(n,nx,ny,dx,dy,v)
	implicit real*8(a-h,o-z)
	dimension v(n),an(100)
	modes=100
	pi=acos(-1.0)
	do 10 in=1,modes
		a1=cos(0.75*in*pi) - cos(0.25*in*pi)
		an(in)=-2.0*a1/(in*pi*sinh(in*pi))
   10	continue
	do 20 j=1,ny
		y=j*dy
		do 18 i=1,nx
			x=dx*i
			s=0
			do 16 jn=1,modes
				s=s+an(jn)*sinh(jn*pi*y)*sin(jn*pi*x)
   16			continue
			ll=(j-1)*nx+i
			v(ll)=s
   18		continue
   20	continue
	return
	end